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Abstract.  A new  class  of  quadrature  rules  for  the  integration  of  both  regular  and  singular  functions 
is  constructed  and  analyzed.  For  each  rule  the  quadrature  weights  are  positive  and  the  class  includes 
rules  of  arbitrarily  high-order  convergence.  The  quadratures  result  from  alterations  to  the  trapezoidal 
rule,  in  which  a small  number  of  nodes  and  weights  at  the  ends  of  the  integration  interval  are  replaced. 

The  new  nodes  and  weights  are  determined  so  that  the  asymptotic  expansion  of  the  resulting  rule, 
provided  by  a generalization  of  the  Euler-Maclaurin  summation  formula,  has  a prescribed  number  of 
vanishing  terms.  The  superior  performance  of  the  rules  is  demonstrated  with  numerical  examples  and 
application  to  several  problems  is  discussed. 
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1.  Introduction 

The  well-known  Euler-Maclaurin  summation  formula  provides  an  asymptotic  expansion  for  the 
trapezoidal  rule  applied  to  regular  functions.  While  the  constant  term  of  the  expansion  is  an  integral, 
the  other  terms  depend  on  the  integrand’s  derivatives  at  the  endpoints  of  the  interval  of  integration. 
This  expansion  is  often  used  to  “correct”  the  trapezoidal  rule  to  a quadrature  with  high-order  conver- 
gence, through  the  use  of  either  known  derivative  values  or  their  finite-difference  approximations. 

In  this  paper  we  derive  asymptotic  expansions,  analogous  to  the  Euler-Maclaurin  formula,  for  func- 
tions with  known  singularities.  In  particular,  functions  with  power  or  logarithmic  singularities  are 
treated.  The  Euler-Maclaurin  formula  and  the  new  asymptotic  expansions  are  used  to  construct  quad- 
rature rules  of  arbitrary  order  convergence.  Each  quadrature  is  constructed  by  altering  the  trapezoidal 
rule:  a few  of  the  nodes  and  weights  at  the  interval  endpoints  are  replaced  with  new  nodes  and  weights 
determined  so  as  to  annihilate  several  terms  in  the  asymptotic  expansion.  The  nodes  always  lie  within 
the  interval  of  integration  and  the  weights  are  always  positive. 

For  a regular  function  / : [0,  1]  -»  K,  we  approximate  /0'  f(x ) dx  with  the  quadrature 

Tn(f)  = h[wi  f(x\h)  + w2  f(x2h ) -4 f wj  f(Xjh ) 

+ f(ah)  + f(ah  + h)  - 1 f /(I  - ah) 

+ wj  f (1  ~ Xjh)  -| f tu  i /(I  - xih)\.  (1) 

There  are  n “internal”  nodes  with  spacing  h = 1 /(n  + 2a  — 1)  and  j “endpoint”  nodes  at  each  end, 
with  the  endpoint  nodes  x\ , . . . , xj  and  weights  w\ , . . . ,Wj  chosen  so  that  the  asymptotic  expansion 
of  T„  as  n — > oo  has  2 j vanishing  terms  and 

fn(f)  = ['  f(x)dx  + 0(h2'+l)  (2) 

Jo 

(Theorem  3.1  and  Corollary  3.2).  The  parameters  a and  j,  and  the  nodes  xi, ...  ,Xj  and  weights 
w i, ...  , Wj,  are  independent  of  n.  The  nodes  and  weights  are  determined  by  2j  nonlinear  equations, 
which  have  a unique  solution,  with 

0 < Xi  < a,  Wi  > 0,  i = , j,  (3) 

provided  a is  sufficiently  large  (Theorem  4.7).  For  integrands  that  are  singular  at  one  endpoint,  Tn 
is  altered  so  that  the  nodes  and  weights  at  that  end  differ  from  those  at  the  other  end  and  depend  on 
the  singularity  (Theorem  3.4  and  Corollary  3.6;  Theorem3.7  and  Corollary  3.8).  For  improper  inte- 
grals in  which  the  integrand  is  oscillatory  and  slowly  decaying,  Tn  is  combined  with  Gauss-Laguerre 
quadrature  to  give  rules  with  high-order  convergence  (Theorem  3.9  and  Corollary  3.10). 

Several  authors  have  studied  the  problems  treated  here.  It  has  been  observed  that  endpoint  correc- 
tions can  be  derived  for  singular  integrands;  Rokhlin  [1]  implemented  such  a scheme  for  integrands 
with  a known  singularity  at  an  interval  endpoint.  He  derived  corrections  to  the  trapezoidal  rule  by 
placing  additional  quadrature  nodes  near  the  endpoint,  with  the  corresponding  weights  determined 
so  that  low-order  polynomials  and  the  singularity  times  low-order  polynomials  were  integrated  ex- 
actly. He  showed  that  under  fairly  general  conditions,  these  weights  had  limiting  values  (up  to  scale) 
as  the  number  of  nodes  in  the  trapezoidal  rule  increased  without  bound  and  that  these  limiting  weights 
could  be  used  to  form  quadrature  rules  with  good  convergence.  Unfortunately,  the  order  of  conver- 
gence of  these  rules  is  restricted  in  practice  by  the  fact  that  the  weights  increase  in  magnitude  rapidly 
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as  the  order  increases.  Efforts  by  Starr  [2]  and  subsequently  Alpert  [3]  reduced  the  growth  in  size  of 
the  weights  with  order,  primarily  by  using  more  weights  than  the  number  of  equations  satisfied  and 
minimizing  their  sum  of  squares.  In  another  approach,  Kress  [4]  uses  all  quadrature  weights  in  the 
quadrature  rule,  rather  than  a few  near  the  endpoints,  to  handle  the  singularity.  More  recently,  Kapur 
and  Rokhlin  [5]  successfully  constructed  rules  of  arbitrary  order  by  separating  the  integrand’s  regular 
and  singular  parts  and  allowing  some  quadrature  nodes  to  lie  outside  the  interval  of  integration. 

The  present  approach  does  not  suffer  from  limitations  on  order  of  convergence,  separation  of  the 
integrand  into  parts,  or  quadrature  nodes  outside  the  interval  of  convergence.  On  the  other  hand,  the 
quadrature  nodes  near  the  interval  endpoints  are  not  equispaced.  Also,  the  equations  for  the  nodes 
x\ , . . . , xj  and  weights  w\, . . . ,Wj,  in  addition  to  being  nonlinear,  are  poorly  conditioned;  the  con- 
ditioning deteriorates  rapidly  with  increasing  order.  Nevertheless,  we  are  able  to  use  an  algorithm 
developed  recently  by  Ma,  Rokhlin,  and  Wandzura  [6]  for  computing  generalized  Gaussian  quadra- 
tures to  obtain  accurate  quadrature  nodes  and  weights.  The  author  would  also  like  to  credit  that  paper 
for  inspiring  the  present  work. 

The  paper  is  organized  around  Section  3,  where  the  new  asymptotic  expansions  are  derived  and 
the  quadratures  defined,  and  Section  4,  where  it  is  shown  that  the  equations  defining  the  quadratures 
actually  have  solutions,  which  are  unique.  These  sections  are  preceded  by  mathematical  preliminaries 
and  followed  by  a discussion  of  the  computation  of  the  quadrature  nodes  and  weights.  Numerical 
examples  are  presented  in  Section  6 and  we  conclude  with  some  applications  and  a summary. 

2.  Mathematical  Preliminaries 

The  material  in  this  section,  which  is  found  in  standard  references,  is  used  in  the  subsequent  devel- 
opment. 

2. 1 . Bernoulli  Polynomials.  The  Bernoulli  polynomials  are  defined  by  generating  function  (see,  for 
example,  Abramowitz  and  Stegun  [7]  23.1.1), 


from  which 


The  Bernoulli  polynomials  satisfy  the  difference  formula 


(4) 


the  differentiation  formula 


B'n(x)  = n Bn_i(x),  n = 1,2,..., 


(5) 


and  the  expansion  formula 


n = 0,  1, . . . . 


(6) 
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2.2.  Euler-Maclaurin  Summation  Formula.  For  a function  / e CP(R)  with  p > 1,  the  Euler- 
Maclaurin  summation  formula  (see,  for  example,  Abramowitz  and  Stegun  [7]  23. 1 .30)  can  be  derived 
by  repeated  integration  by  parts.  We  first  consider  the  interval  [c,  c + h]  and  apply  (5)  to  obtain 

-c+h 


/c+h  c 1 

f(x)dx=h  j Bq(\  — x)  f(c  + xh)  dx 


= -h 


Bx(\-x) 

1! 


f(c+xh ) 


+ h 


f 


Bi(\-x) 

1! 


f\c  +xh)dx 


= -E 


hr+l  Br+i(\  -x) 

(r  + 1)! 


f(r\c  + xh) 


1 


0 


where  we  have  used 


+ h~‘  f'w-*' 

Jo 


P'- 


/(/>)(c  + xh ) dx 


= * _ g [/M(e + » - /«(c)] 


+ *<’+'  ['  f^(c  + xh)dx, 

Jo  p'- 


-B,(0)  = 51(l)  = i 

B„(0)  = fi„(l)  = 5n,  n#l. 


(7) 


To  derive  the  Euler-Maclaurin  formula  for  the  interval  [a,  £>],  we  let  h = (b  — a)/ n and  c = a + ih 
in  (7),  sum  over  / = 0,  1 , . . . , n — 1 , and  rearrange  terms,  to  obtain 


f(a)  , */  , , , , fib)' 

— b f(a+h)-\ h fib  - h)  + — — 


= Ja  f^dx  + 5Z  V+T)^  [/(r)(Z?)  “ /(r,(a)3 

-^+1  f /*>(?  + ih+xh) 

Jo  P'-  l to 


rf*.  (8) 


The  expression  on  the  left-hand  side  of  (8)  is  the  well-known  trapezoidal  rule.  Evaluation  of  the  ex- 
pression on  the  right-hand  side  of  (8)  is  simplified  by  the  fact  that  Z?2r+i  = 0 for  r > 1. 


2.3.  Generalized  Riemann  Zeta-Function.  The  generalized  Riemann  £ -function  is  defined  by  the 
formula 


S(s,  v) 


OO 


E 


1 

( v + n)s  ' 


Re(s)  >1,  r^O,  — 1, . . . . 


This  function  has  a continuation  that  is  analytic  in  the  entire  complex  5 -plane,  with  the  exception  of 
s = 1 , where  it  has  a simple  pole.  In  what  follows,  we  shall  be  concerned  primarily  with  real  s and 
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v,  with  s < 1 and  v > 0.  We  will  use  the  following  representation  derived  from  Plana’s  summation 
formula  (see,  for  example,  Erdelyi  [8]  1.10  (7)): 


v ) = 


, 1—  S 


S — 1 


+ 


-+2v  i 


sin(j  arctant)  dt 
(1  + t2y/2  e2nvt  - 1 


Re(u)  > 0. 


(9) 


Equation  (9)  can  be  used  to  derive  the  asymptotic  expansion  of  £ as  v -*  oo.  We  treat  the  integral 
as  a sum  of  Laplace  integrals,  each  with  an  asymptotic  expansion  given  by  Watson’s  lemma  (see,  for 
example.  Bender  and  Orszag  [9]  p.  263),  and  obtain 

,l-.r  v~s  1 P 


Vl~S  V~S  1 ^ 

f(j,„)  = — + 1-  + — £ 


s + 2r  — 2 
2 r 


Blr 


V 


s+2r—  1 


+ 0(v 


—s—2p— 


(10) 


as  v oo,  with  s € C,  s 1,  and  p an  arbitrary  positive  integer.  Equation  (10)  is  a slight  general- 
ization of  [8]  1.18  (9).  There  is  a direct  connection  between  the  Bernoulli  polynomials  and  f , 


Bn(x) 


= -£(1  -n,x),  n = 1,2, ...  , 


and  generalizations  of  the  difference  and  differentiation  formulae  hold: 

it- 1 , 


£(s,  v)  - Us,  V + k)  = V - — — — , 

jz o (.v  + iy 

d£(s,  v ) 


dv 


= -s  Z(s  + 1,  v). 


(11) 

(12) 


2.4.  Orthogonal  Polynomials  and  Gaussian  Quadrature.  Suppose  that  co  is  a positive  continuous 
function  on  the  interval  (a,  b ) and  co  is  integrable  on  [a,  b].  We  define  the  inner  product  with  respect 
to  co  of  real- valued  functions  / and  g by  the  integral 

f(x)g(x)co(x)dx. 

J a 

There  exist  polynomials  po,  p\,  P2, . . . , of  degree  0,  1,2,...,  respectively,  such  that  ( pn , pm ) = 0 
for  n ^ m (orthogonality);  they  are  unique  up  to  the  choice  of  leading  coefficients.  With  leading  co- 
efficients one,  they  can  be  obtained  recursively  by  the  formulae  (see,  for  example,  Stoer  and  Bulirsch 
[10],  p.  143) 

Po(x)  = 1 (13) 

Pn+lix)  = (X  -Sn+\)pn(x)  ~ yn+]2  pn-i(x),  U > 0,  (14) 


where  p-\  (x)  = 0 and  8n,  yn  are  defined  by  the  formulae 

<Wl  = (xpn,  Pn)/ (Pn  i Pn),  n > 0, 


Yn+l2 


fo  n = 0, 

\(Pn,  Pn)/(Pn-l,  Pn- 1)  n>\. 


The  zeros  x”, . . . , x”  of  pn  are  distinct  and  lie  in  the  interval  ( a , b).  There  exist  positive  numbers 
, . . . , oonn  such  that 


f(x)co(x)dx  = ^00/  f(x/) 

i=i 


(15) 
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whenever  / is  a polynomial  of  degree  less  than  2 n.  These  Christoff  el  numbers  are  given  by  the  for- 
mula (see,  for  example,  Szego  [11]  p.  48) 


(Pn— 1 * Pn—  1 ) 

Pn-iW)  p'n(x?)’ 


(16) 


Moreover,  if  (o{x ) = (b  — x)  r(x)  with  r integrable  on  [ a , b],  then,  with  the  definition  *"+1  = b, 
there  exist  positive  numbers  r", . . . , r"+]  such  that 


(17) 


whenever  / is  a polynomial  of  degree  less  than  or  equal  to  2 n.  These  modified  Christoffel  numbers 
are  given  by  the  formula 


(O' 


l 

Pn(x) 
Pn  iff) 


r(x)  dx 


i = 1, . ..  ,n, 
i = n + 1 , 


(18) 


where  oJ]  is  given  by  (16). 

The  summation  in  (15)  is  the  n node  Gaussian  quadrature  with  respect  to  a>,  while  that  in  (17)  is  an 
n + 1 node  Gauss-Radau  quadrature  with  respect  to  r. 


3.  Hybrid  Gauss-Trapezoidal  Quadrature  Rules 

In  this  section  we  introduce  new  quadrature  rules  for  regular  integrands,  singular  integrands  with  a 
power  or  logarithmic  singularity,  and  improper  integrals,  and  determine  their  rate  of  convergence  as 
the  number  of  quadrature  nodes  increases. 

For  notational  convenience  we  generally  consider  quadratures  on  canonical  intervals,  primarily 
[0,  1].  It  is  understood  that  these  are  readily  transformed  to  quadratures  on  any  finite  interval  [a,  b ] 
by  the  appropriate  linear  transformation  of  the  nodes  and  weights. 

3.1.  Regular  Integrands.  For  j,  n positive  integers  and  a € R+  = {*  € R|;c  > 0},  we  define 
a linear  operator  Tff  on  C([0,  1]),  depending  on  nodes  x\,  ...  , x}  and  weights  wi, ...  , Wj,  by  the 
formula 

j n- 1 j 

Tff(f ) = h^Wi  f (Xjh)  +h^2  f(ah  + ih)  +h^2wi  /(I  - xth),  (19) 

i=l  i=0  i=l 

where  h = (n  + 2a  — l)-1  is  chosen  so  that  ah  + {n  — \)h  = 1 — ah. 

Theorem  3.1.  Suppose  f e Cp([0,  1]).  The  asymptotic  expansion  ofTff(f)  as  n -*  oo  is  given  by 
the  formula 


Tff(f) 


= /' 


f (x)  dx 


p- 1 


+ Y.h 


-]Y  fPh 


r= 0 


r\ 


(20) 


r + 1 


HYBRID  GAUSS-TRAPEZOIDAL  QUADRATURE  RULES 


7 


Proof.  We  apply  the  Euler-Maclaurin  formula  (8)  on  the  interval  [ah,  1 — ah]  to  obtain 


h Y\  f(ah  + ih)  = h fiah)  + ^(1 — — + [ f(x)  dx 
7=o  2 Jah 


f(ah)  + f(\  - ah) 


• l— ah 


('  + D! 

We  now  combine  (19)  and  (21),  the  equality 


P~  * Lr+^  J) 

+ T.  iVTTTT  h<r)(1  ~a  h)-  fr\ah)]  + 0(h”+').  (21) 


p\—ah  p 1 pah  p I 

/ f(x)dx=  / / (x)dx  — / / (x)dx  — / f(x)dx, 

Jah  Jo  Jo  J 1 —ah 

Taylor  expansion  of  all  quantities  about  h = 0,  the  Bernoulli  polynomial  expansion  formula  (6),  and 
difference  formula  (4)  to  obtain  (20).  □ 


Corollary  3.2.  Suppose  the  nodes  jcj  , . . . ,xj  and  weights  w\ , . . . ,Wj  satisfy  the  equations 

r Br+l(a)  0-1 

> WiXi  = — , r = 0,  1, . . . , 2]  - 1. 

r + 1 


(22) 


i=i 


Then  Tfa  is  a quadrature  rule  with  convergence  of  order  2 j + 1 for  f € Cp([  0,  1])  with  p > 2 j. 
Moreover, 


(23) 


as  n — » oo,  provided  /(2-,)( 0)  4-  /(2v)(l)  ^ 0. 


Corollary  3.3.  Suppose  Xj  = a — 1 and  the  remaining  nodes  x\,  . . . , Xj-\  and  weights  w . ,Wj 
satisfy  the  equations 


TwiXjr  = 


i=i 


Br+\{a) 

r + 1 


r =0,1,...  ,2; -2. 


(24) 


Then  Tfa  is  a quadrature  rule  with  convergence  of  order  2 j for  f € Cp([0,  1])  with  p > 2 j — 1 . 
Moreover, 

**</)  - j[‘  /(*)*  ~*^ j - M4)  (25, 

as  n ^ oo,  provided  /(2-'~1)( 0)  — /(2j,_1)(l)  # 0. 

We  shall  see  below  that  (22)  has  a solution  with  the  nodes  and  weights  all  positive  if  a is  sufficiently 
large  and  that  numerical  solution  of  (22)  is  equivalent  to  computing  the  roots  of  a particular  polyno- 
mial. This  statement  holds  for  (24)  as  well. 


8 


BRADLEY  K.  ALPERT 


3.2.  Singular  Integrands.  For  j,  k,  n positive  integers  and  a,b  g R+,  we  define  a linear  operator 
S}nkab  on  C((0,  1]),  depending  on  nodes  Vi, . . . ,Vj,x  i, ...  ,xk  and  weights  u\, . . . , Uj,  w\, . . . , wk, 
by  the  formula 


j n—  1 k 

SJnkab(g ) -h^Ui  g(Vih)  + h^  S(ah  + ih)  + h ^ w(  g(l  - xth),  (26) 

i'=l  i=0  i=l 

where  h = (n  + a + b — l)-1  is  chosen  so  that  ah  + (n  — \)h  = 1 — bh. 

Theorem  3.4.  Suppose  g(x)  =xyf(x),  where  y > —1  and  f e Cp([0,  1]).  The  asymptotic  expan- 
sion ofSJnkab(g)  as  n — » oo  is  given  by  the  formula 


r * p~\  f(n( n\  r j 

SJnkab(g)=  / g(x)dx  + ^hy+r+l j — Viy+r  + f(-y  - r,  a) 

Jo  r=o  r ■ ‘ i=l 

+ g^(-1)rf)(1)  ’ 


r= 0 


Y w,  x,'  - + O(/ip+1+minf0'J'1)-  (27) 

sr  r + 1 I 


Proof.  For  c € M+,  we  define  polynomials  pcQ , pf, . . . , in  analogy  with  the  Bernoulli  polynomials, 
by  the  formula 


>:;<*>  = E 


r=0 


C(-y-r,c)(l-c-xr-r. 


(28) 


Differentiating,  we  verify  that 

—pn(x)  = -n  Pn-fx),  n = 1,2,... 
and  combining  the  £ difference  formula  (11)  with  (28)  we  obtain 


pS(1)-^+'<0)  = 


0 


72  = 0, 
72  > 1. 


Additionally,  we  define  functions  qc0,  q\, . . . by  the  formula 


$£(*)  = 


(x  + c)^ 

(-!)"(*  +c)y+n 
(y  + l)  • • • (y  + 72) 


Pn-\(X) 

(n-  1)! 


72  = 0, 

#1  = 1,2,..., 


(29) 


and  observe  that 


-^qcn(x)  = -qcn_  fx), 


72  = 1,  2,  ...  . 


HYBRID  GAUSS-TRAPEZOIDAL  QUADRATURE  RULES 


9 


With  these  definitions,  the  proof  follows  that  of  the  Euler-Maclaurin  formula. 


r 1-bh  n-2 

/ xyf(x)dx 

J ah 


q0  (x)  f(ah  + ih  + xh ) dx 


n-2  r\ 

= Z»'+r 

i=0  •'0 

n-2  i p- 1 I 

= E - E"V+r+1C+;U)  /<r)(aA  + + *A) 

1=0  l r=0  o 

_j_  /zy+p+1  t q“+‘(x)  f{p)(ah  +ih  + xh)dx  j 

n-2 

= h ^^{ah  + ih)Y  f(ah  + z/z) 

i=i 
p- 1 

- [9"rJ(l)/w(l  - bh)  - 9r“+,(0)/w (<,*)] 


+ *y+p+| 


/•l  f n-2 

X Is 


<Jp+'(x)  f{p\ah  + ih  + xh)  \dx. 


(30) 


Taylor  expansion  of  f(r){ah)  about  h = 0,  the  definitions  (28)  and  (29)  for  pcn  and  qcn,  and  the  binomial 
theorem  combine  to  yield 


p-\ 

£Ay+'+y+1(0)/w(a/0  = 


r=0 


Pz\  hy+k+lfw(0) 

k= 0 


k\ 


aY+k+\ 


_L_  b-  _1_  1 


+ C(~K  -k,a  + 1) 


+ 0(hY+p+1).  (31) 


Likewise,  Taylor  expansion  of  /(r)(l  — bh)  about  h = 0,  the  definitions  for  pcn  and  qcn,  the  asymp- 
totic expansion  (10)  for  £,  the  Bernoulli  polynomial  expansion  formula  (6),  and  the  binomial  theorem 
combine  to  yield 


p-\ 


]TV+r+1*2?-2(l)  /«( l - bh)  = 


r= 0 


g^_+1(-i)^w(i) 


*=0 


k! 


bk+1  Bk+l(b+  1) 
* + 1 Jfc  + 1 


+ 0(hp+1).  (32) 


We  now  combine  (26)  and  (30)  through  (32),  the  equality 


rl  —bh  rl  rah  r 1 

/ g(x)<ix  = / g(x)dx-  / g(x)dx-  / g(x)dx, 
J ah  JO  Jo  J \—bh 


expansion  of  the  latter  two  integrals  about  h = 0,  and  the  difference  formula  (4)  for  the  Bernoulli 
polynomials  and  ( 1 1 ) for  £ to  obtain  (27).  □ 
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Corollary  3.5.  Suppose  the  nodes  u\, ...  , uj  and  the  weights  v\,  . . . , Vj  satisfy  the  equations 


(33) 


and  the  nodes  x\ , . . . , Xj  and  weights  w\ , . . . , wj  satisfy  the  equations 


r = 0,  1 , . . . , 2j  — 1 . 


(34) 


Then  Sjnjab  is  a quadrature  rule  with  convergence  of  order  2j  + \+  min{0,  y } for  g,  where  g(x)  = 


xy  f(x),  with  f € Cp([ 0,  1])  and  p >2j.  Moreover, 


(35) 


Y >0,  g(2j\  1)#0, 


as  n — > oo. 


Corollary  3.6.  Suppose  the  nodes  u i , . . . ,Uj  and  the  weights  v\ , . . . ,Vj  satisfy  the  equations 


j 


]pH«  v/+r  = -£(-y  — r,  a),  r = 0,  1, . . . , j - 1, 


i=i 


(36) 


and  the  nodes  x\ , . . . ,Xk  and  weights  w\, . . . , Wk  satisfy  the  equations 


r = 0,  1,...  ,2k  - 1. 


(37) 


Then  S]nkah  is  a quadrature  rule  with  convergence  of  order  min  {7  + 1,  y + j + 1 , 2k  + 1 } for  g , where 


g(x)  = xyf{x)  + ffx),  with  <p,ijs  e Cp([0,  1])  and  p > min{y,  2k}. 

In  Corollaries  3.5  and  3.6  an  even  number  of  constraints  on  the  nodes  and  weights  at  both  ends  of  the 
interval  are  considered.  Clearly,  there  are  analogous  quadrature  rules  arising  from  an  odd  number  of 
constraints  at  one  or  both  ends;  these  are  similar  and  explicit  presentation  of  them  is  omitted.  We  now 
consider  a different  type  of  singularity. 
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Theorem  3.7.  Suppose  g(x)  = f (x)  log x,  where  f e Cp([0,  1]).  The  asymptotic  expansion  of 
SJnkab(g)  as  n — > oo  is  given  by  the  formula 


SJ„tab(g)  = f 

Jo 


g(x)  dx 


, ^,r+i/(r)(0)  [ r i / / \ yi,  \ Br+\(a.) 

+ j — j 2^ui  vi  log (Vjh)  - C (~r,a) ——log h 

r= o r ■ y i=\  r + l 


p-\ 


, r Br+l(b)\  +1 

+ 2^ | 2^  ~ r+1  | + Q(^P+1l0g/2), 


(38) 


r=0 


where  f ' denotes  the  derivative  oft;  with  respect  to  its  first  argument. 


Proof  This  asymptotic  expansion  is  derived  from  that  of  Theorem  3.4  by  differentiating  (27)  with 
respect  to  y and  evaluating  the  result  at  y — 0.  □ 


Corollary  3.8.  Suppose  the  nodes  u\, . . . , u t and  the  weights  Uj, . . . , Vj  satisfy  the  equations 


J 

Y2  ui  vir  log  v‘  = f'C-r.a),  r = 0, 1,...  , j - 1, 


i=i 


r Br+\(a) 

> UjVi  = — — , r = 0,  1, . . . , J — 1, 

U r + 1 

and  the  nodes  X\ , . . . ,xk  and  weights  w\ , . . . , satisfy  the  equations 


Y2  Wi  xf  = 


Br+l(b) 

r + 1 


r = 0,  1 , . . . , 2k  — 1 . 


(39) 


(40) 


Then  SJnkab  is  a quadrature  rule  with  error  of  order  0(rmn{hj+l  log  h,  h2k+l})for  g,  where  g(x)  = 
4>(x)  log*  + rj/(x),  with  0,  \f  € Cp([ 0,  1])  and  p > min{y,  2k}. 

3.3.  Improper  Integrals.  For  j,  n positive  integers,  we  define  a linear  operator  RJn  on  C([n,  oo)), 
depending  on  nodes  Jti , . . . , xj  and  weights  wi , . . . , wj,  by  the  formula 


j 

Rjn(g)  = J2wkS(n  +xk). 


k=  1 


(41) 


Theorem  3.9.  Suppose  that  g(x)  = eiyx  f(x),  where  y € R,  y ^ 0,  and  f € Cp{[  1,  oo)),  and  that 
there  exist  positive  constants  f,ao,...  , ap,  such  that 

|/(r)(*)|  < -^77  forx>  1,  r = 0, . . . , p.  (42) 

The  asymptotic  expansion  of  RJn(g)  as  n — - oo  is  given  by  the  formula 


K(g) 


-r 


g(x)  dx 


+ eiyn 


E 


r= 0 


f{r)(n) 

r\ 


Y2wkxkreiyXk 
k=  1 


r+1 


+ 0(n~p).  (43) 
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Proof.  We  integrate  by  parts  repeatedly  to  obtain 


p- 1 


/oo  y 1 . : . r+j  />oo  /i\p 

e,rxf(x)dx=e‘'"‘J2{-)  /W(")+j  /wW*, 

and  in  (41)  we  compute  the  Taylor  expansion  of  / about  x*  = 0 to  get 

P_1  /(r)(")  r , f^in+Zik) 


j—  r f^Un) 

Rite)  = Ec"''*u’<  E ;L- + 

*=1  1 r=0  r • 


(44) 


(45) 


where  £*  lies  between  0 and  xk  for  k = l, , j.  Now  combining  (42),  (44),  and  (45)  we  obtain 
(43).  □ 


Example.  The  function  / defined  by  the  formula 


00 


/w  = E 

r=0 


ar 

xP+r  ’ 


(46) 


with  ^ |or|  < oo,  satisfies  the  assumptions  of  Theorem  3.9  for  every  positive  integer  p.  We  remark 
that  Theorem  3.9  can,  in  some  instances,  be  generalized  to  y = 0,  but  the  corresponding  asymptotic 
expansion  depends  on  a more  detailed  knowledge  of  /.  For  / given  by  (46),  for  example,  the  quad- 
rature nodes  and  weights  for  y = 0 depend  on  f. 


Corollary  3.10.  Suppose  f e Cp{[  1,  oo)),  / satisfies  (42)  for  x € R,  and  f is  analytic  in  the  half- 
plane Re(;c)  > a for  some  a € M.  Suppose  further  that  iq, . . . ,Vj  are  the  roots  of  the  Laguerre 
polynomial  Lj  of  degree  j,  that  coefficients  u\, ...  ,uj  satisfy  the  equations 

j 

y ^ ukvkr  e~Vk  — r!,  r = 0,  1, 1,  (47) 

*=i 

and  that  the  operator  RJn  is  defined  with  nodes  xk  = (/ /y)vk  and  weights  wk  = (i /y)uk  for  k = 
1,  . . . , j . Suppose  finally  that  Tm  is  defined  to  be  the  quadrature  rule  with  nodes  and  weights 
satisfying  (22),  but  translated  and  scaled  to  the  interval  [1,  n].  Then  for  p > 2 j,  the  expression 
T(n-\)n(g)  + Rn(g)  « tin  approximation  for  the  integral  g(x)dx,  where  g(x)  = e‘YX  f(x),  with 
error  of  order  0(n~ 2i)  as  n — > oo. 


Proof  This  result  is  just  a combination  of  the  quadrature  rule  of  Corollary  3.2,  for  the  interval  [1 , n], 
with  the  asymptotic  expansion  of  Theorem  3.9,  for  the  interval  [n,  oo),  provided 


Yukvk'e  Vk  = r !, 

k= 1 

But  (48)  follows  from  (47),  the  equations 

POO 

r\=  xre~x 
Jo 


dx. 


r = 0,  1, . . . , 2j  — 1. 


r = 0,  1, . 


(48) 


f 


Lj(x)  Lk(x)  e x dx  = 0 for  j k , 


(see,  for  example,  [7]  6.1.1  and  22.2.13)  and  the  fact  that  Gaussian  quadratures,  which  are  exact  for 
polynomials  of  degree  less  than  twice  the  number  of  nodes,  have  nodes  that  coincide  with  the  roots 
of  the  corresponding  orthogonal  polynomials  (see  Section  2.4).  □ 
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We  have  completed  the  definition  of  the  new  quadratures,  along  with  the  demonstration  of  their  as- 
ymptotic performance.  We  shall  see  that  the  existence  of  these  rules,  which  depends  on  the  solvability 
of  the  nonlinear  systems  of  equations  that  define  the  nodes  and  weights,  is  assured  by  the  theory  of 
Chebyshev  systems.  The  uniqueness  of  the  rules  is  similarly  assured.  These  issues  of  existence  and 
uniqueness  are  treated  next. 


4.  Existence  and  Uniqueness 


4. 1 . Chebyshev  Systems.  Material  of  this  subsection  is  taken,  with  minor  changes,  from  Karlin  and 
Studden  [12].  Suppose  / is  an  interval  of  K,  possibly  infinite.  A collection  of  n real-valued  continuous 
functions  fi,  ...  , fn  defined  on  I is  a Chebyshev  system  if  any  linear  combination 

n 

fix)  = ^Qifiix), 

i=l 


with  a.i  not  all  zero,  has  at  most  n — 1 zeros  on  /.  This  condition  is  equivalent  to  the  statement  that 
for  distinct  x\ , . . . , xn  in  I, 


det 


( f\{X\)  •••  fnix  i)> 
\f\iXn)  fniXn)) 


#0. 


(49) 


The  Chebyshev  property  is  a characteristic  of  the  space,  rather  than  the  basis:  if  f\ , . . . , /„  is  a Cheby- 
shev system,  then  so  is  any  other  basis  of  span{/i , . . . , /„ }.  If  u is  a continuous,  positive  function  on 
/ then  scaling  by  u preserves  a Chebyshev  system.  Finally,  if  u is  strictly  increasing  and  continuous 
on  interval  J with  range  /,  then  f\ou, ...  , fn  ou  is  a Chebyshev  system  on  J if  and  only  if  fi, ...  , fn 
is  on  I.  (Here  /,  o u denotes  the  composition  u followed  by  fi.) 

The  best-known  example  of  a Chebyshev  system  is  the  set  of  polynomials 


l,x, ...  ,xn  1 

on  any  interval  / cR.  We  shall  be  concerned  also  with  the  Chebyshev  systems 

l Xy  X Xy  + l ^(n-U/2  xy+(n- 1)/2 

1,  log*,  x,  x log*, . . . , x{n~1)/2,  x(n~1)/2  log* 

on  / = (0,  a],  where  y € K\Z  and  a > 0.  These  systems  are  special  cases  of  the  system  of  Muntz 
functions  (see,  for  example,  Borwein  and  Erdelyi  [13]  p.  133) 

M = { xy‘  log*  x | k = 0, . . . , ni  — 1 , i = 1 , . . . , j } (50) 


on  I = (0,  oo),  where  yi, ...  ,yj  are  distinct  real  numbers  and  n\, ...  , nj  are  positive  integers  with 
= n.  To  see  this  is  a Chebyshev  system,  suppose  / g span  M and  use  induction  in  n on 
(d / d log x)[ f(x)  x~yj],  in  combination  with  Rolle’s  theorem.  Another  Chebyshev  system  that  will 
arise  is  the  system 

L = {(*  + yi)~k  | k = 1, . . . ,nit  i = 1, . . . , j } (51) 


on  / = [0,  oo),  where  yi, . . . , yj  are  distinct  positive  real  numbers  and  n\, . . . ,n}  are  positive  in- 
tegers with  Jfni  = n.  This  is  indeed  a Chebyshev  system,  for  if  / € span  L then  the  function 
fix)  nu  i ix  + Yi)n'  is  a polynomial  in  x of  degree  n — 1. 
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Suppose  /i , . . . , /„  is  a Chebyshev  system  on  the  interval  7.  The  moment  space  Mn  with  respect 
tofu...  , fn  is  the  set 


Mn  = {c  = (ci,  Ci  = j fi(x)da(x),  i = 1 n j, 


(52) 


where  the  measure  o ranges  over  the  set  of  nondecreasing  right-continuous  functions  of  bounded  vari- 
ation on  7.  It  can  be  shown  that  /Mn  is  the  convex  cone  associated  with  points  in  the  curve  C,  where 

C = {(/i(x), ...  , /„(*)) | * € /}. 

In  other  words,  Mn  can  be  represented  as 

f P 1 

= |c  = T>;  yj  olj  > 0 ,yj  eC,  j = 1, . . . , p,  p > 1 [. 

7=i 


The  index  1 (c)  of  a point  c of  Mn  is  the  minimum  number  of  points  of  C that  can  be  used  in  the 
representation  of  c,  under  the  convention  that  a point  {f\  (x), . . . , fn(x ))  is  counted  as  a half  point  if 
x is  from  the  boundary  of  7 and  receives  a full  count  otherwise.  The  index  of  a quadrature  involving 
*i, . . . , xp  is  determined  by  counting  likewise. 

Proofs  of  the  next  three  theorems  are  somewhat  elaborate  and  are  omitted  here;  they  can  be  found 
in  Karlin  and  Studden  [12]. 


Theorem  4,1  ([12]  p.  42).  Suppose  I = [ a,  b ] is  a closed  interval.  A point  c € /Mn,  c ^ 0,  is  a 
boundary  point  of  SMn  if  and  only  if  1(c)  < n/2.  Moreover,  if  cr  is  a measure  corresponding  to  a 
boundary  point  c € !Mn,  then  there  is  a unique  quadrature 

p 

IV i fr{Xi)  = 

i=l 

where  p < (n  + l)/2,  a < x\  < X2  < ■ • ■ < xp  < b,  and  Wi  > 0,  i = 1,  . . . , p. 


fr(x)  da(x). 


r = 1, 


n, 


(53) 


Theorem  4.2  ([12]  p.  47).  Suppose  7 = [a,  b]  is  a closed  interval.  Any  point  c in  the  interior  of  SMn 

satisfies  1(c)  = n /2.  Moreover,  if  a is  a measure  corresponding  to  c,  then  there  are  exactly  two 

quadratures 

Wi  fr(xi)  = J^fr(x)da(x),  r = \,...,n,  (54) 

of  index  n/2,  where  Wi  > 0,  i = 1, . . . , p.  In  particular,  if  n = 2m,  then  p = m or  p = m + \ and 

a < x\  < X2  <■■■<  xm  < b or  (55) 

a — x i < x2  < ■ ■ ■ < xm+\  = b\  (56) 

if  n = 2m  + 1,  then  p = m + \ and 

a = x i < X2  < ■ ■ • < xm+\  < b or  (57) 

a < x i < X2  < ■ ■ ■ < xm+i  = b.  (58) 


E 


Theorem  4.3  ([12]  p.  65).  Let  Tn  denote  the  nonnegative  linear  combinations  of  f\, ...  , fn, 

n 

"Pn  = {/I  fix)  = ^2, ai  ft  (*)  and  /(*)  >0/»rfl/Zjce7}. 

i=i 


(59) 
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The  point  c = (cj , . . . , cn)  is  an  element  ofMn  if  and  only  if 

n n 

yt  at  fi  € Tn  implies  ^ a,  c,  > 0.  (60) 

i=i  1=1 

Moreover,  c is  in  the  interior  of  Mn  if  and  only  if 

n n n 

Tai  fi  G Tn  and  > 0 imply  >0.  (61) 

i=i  i=i  !=i 


Theorem  4.4  ([12]  p.  106).  Suppose  f(x)  = xl  1 for  i = 1,...  ,n,andl  = [a,  b].  If  n = 2m, 
then  c = (ci, . . . , cn)  is  an  element  of  Mn  if  and  only  if  the  two  quadratic  forms 

m m 

y [c,+ j -aci+j- 1 ] a?,- a j and  ^ [bci+J-i  - ci+J ] (62) 

<u=i  <o=i 

are  nonnegative  definite.  Ifn  — 2m  + 1,  then  c € SMn  if  and  only  if  the  two  quadratic  forms 

m+ 1 m 

y Ci+j -idiot j and  £[<*  + b)ci+J  -ab ci+j- 1 - Ci+j+i]  ptpj  (63) 

<o=i  <o=i 

are  nonnegative  definite.  Moreover,  for  either  parity  ofn,  c is  in  the  interior  of  SMn  if  and  only  if  the 
corresponding  quadratic  forms  are  both  positive  definite. 


Proof.  A theorem  of  Lukacs  (see,  for  example,  Szego  [11]  p.  4)  states  that  a polynomial  / of  degree 
n — 1 that  is  nonnegative  on  [a,  b]  can  be  represented  in  the  form 

'(x  -a)p(x)2  + (b  -x)q(x)2  n = 2m, 
p(x)2  + (b  — x)(x  — a)q(x)2  n = 2m  + 1, 


/(*)  = 


where  p and  q are  polynomials  such  that  the  degree  of  each  term  in  (64)  does  not  exceed  n — 1.  The 
combination  of  (64)  and  Theorem  4.3  proves  the  theorem.  □ 

4.2.  Muntz  System  Quadratures.  The  systems  of  equations  (22),  (24),  (36),  and  (39)  that  define  the 
quadrature  rules  of  Section  3 are  special  cases  of  the  system  of  equations 


r«/2i 


y Wmxmy‘  log*  xm  = (-1  )k+lS(k)(-yi,a),  k = 0, . . . , n,  - 1,  i = 1, ...  , j, 


(65) 


m- 1 


for  distinct  real  numbers  yi, ...  , Yj  and  positive  integers  ni, ...  ,nj  with  ni  = «.  Here  £(<:)  de- 
notes the  kth  derivative  of  £ with  respect  to  its  first  argument.  The  existence  and  uniqueness  of  the 
solution  of  (65)  follow  from  the  existence  and  uniqueness  of  quadratures  for  Chebyshev  systems,  once 
it  is  established  that  there  is  a measure  oa  with 

xy‘  log** daa(x)  = (-l)*+1^(*)(-y,,  a),  = 0 m - 1,  i = 1, . . . , j,  (66) 

in  other  words,  that  the  moment  space  IMn  of  the  Chebyshev  system  of  Muntz  functions 

M = { xyi  log* x | k = 0, .. . , n{  - 1,  i = 1, . . . , j),  (67) 

on  (0,  a]  contains  the  point 

C = ( - c(-n,  a), . . . , (-1)»;? (»;-U (-Yj,  a)). 


(68) 
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We  will  show  that  this  condition  is  satisfied  provided  that  a is  sufficiently  large.  It  would  be  convenient 
to  have  tight  bounds  for  a,  in  particular  for  systems  (22),  (24),  (36),  and  (39),  but  it  appears  that  such 
bounds  are  difficult  to  obtain.  Even  for  the  regular  cases  (22)  and  (24),  where  by  Theorem  4.4  the 
existence  of  oa  is  equivalent  to  the  positive  definiteness  of  two  matrices,  precise  bounds  for  arbitrary 
j appear  difficult.  (Numerical  examples  below  provide  evidence  that  a/j  may  be  chosen  as  small  as 
5/6.) 

Theorem  4.5.  Suppose  y\,  . . . ,Yj  are  distinct  real  numbers,  each  greater  than  —1,  and  n i, . . . , n} 
are  positive  integers  with  Y ni  — n-  Then  for  sufficiently  large  a there  exists  a measure  oa  such  that 
the  system  of  equations  (66)  is  satisfied  and  c defined  by  (68)  is  in  the  interior  of  the  moment  space 

Mn. 

Proof.  We  construct  a continuous  weight  function  o'a  satisfying  (66)  and  show  that  for  sufficiently 
large  a , crffx)  is  positive  for  x e [0,  a\. 

We  linearly  combine  the  equations  of  (66)  to  obtain  the  equivalent  system 


J ( x/a)yi  lo  £(x/a)doa(x)  = (-1)*+1  ^ ^ - — "'j"'""'  log*  * a' 


k = 0, ...  ,nt  — 1,  i = 1, .. . , j,  (69) 

where  we  have  used  the  binomial  theorem  to  expand  log k(x/a)  = (log*  — loga)*.  We  define  the 
weight  o’a  by  the  formula 


n—\ 

&a(x)  = (x/af 

m= 0 


and  combine  (69),  (70),  and  the  equalities 


f xy  log* 
Jo 


x dx  = 


(-1  )kk\ 
(l+y)k+r 

to  obtain  the  equations  in  aoM, . . . , 


y > -1,  k = 0, 1, 


tS(l  +Y,+mY*'  ’ r)  «'+»  g 


r= 0 


(70) 


k =0, ...  ,nt  - 1,  i = 1, ....  7.  (71) 


This  n -dimensional  linear  system  is  nonsingular,  since  the  set  of  functions 

{(a  + Yi  + 1)-*|  k = 1, . . . , «j,  / = 1, ....  7} 

forms  a Chebyshev  system  on  [0,  00),  as  established  at  (51).  Thus  (71)  possesses  a unique  solution 

&0 .ai  • • • 1 Otn-l.a- 

We  now  determine 


oim  — lint  cim  a , 

cl — >00 


m = 0,  . . . n — 1. 


(72) 
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The  asymptotic  expansion  of  £(r)(— y,-,  a)  as  a — ► oo  can  be  derived  by  differentiating  (10);  the  first 
several  terms  are  given  by 


C<r>(-y;,a)  = a'+KT 


^ /r  \ (—  1 ),+1  (r  — /) ! log^  a 


1=0 


(1  + Yi) 


\r— /-f 1 


+ log7-  a). 


(73) 


Combining  (7 1 ) and  (73),  changing  the  order  of  summation,  and  twice  applying  the  product  differen- 
tiation rule 

^7 (f(Y)g(Y))  = £ Q f{s)(Y)8ir~s)(Y) 


we  obtain 


am(-\)kk\  _ (~\)kk\ 

A(1  +Yi+*n)k+l  0 +Yi)k+r 


k = 0,  . . . , Hi  — 1 , 7 = 1, 


7, 


which  immediately  reduces  to 


a„,  = 


| 1 m — 0, 

I 0 m = 1, . . . , n — 1. 


(74) 


The  combination  of  (70),  (72),  and  (74)  gives 

lim  o'  (ax)  = 1,  x e [0,  1], 

<2— XX) 

which  implies  that  for  a sufficiently  large,  o'a  (x)  > 0 for  a € [0,  a].  The  point  c defined  by  (68)  is  in 
the  interior  of  !Mn,  since  small  perturbations  of  c will  preserve  the  positivity  of  o'a.  □ 


Theorem  4.2  of  Section  4.1  ensures  the  existence  of  Gaussian  quadratures  for  a Chebyshev  system 
/i,...  , fn  defined  on  an  interval  /,  under  the  assumption  that  I is  closed,  whereas  the  system  M of 
(67)  is  Chebyshev  on  / = (0,  a].  As  a consequence,  we  require  the  following  result. 

Theorem  4.6.  Suppose  the  collection  of  functions  f\ , . . . , fn  forms  Chebyshev  system  on  I = (a,  b\ 
and  each  is  integrable  on  [ a , b]  with  respect  to  a measure  o corresponding  to  a point  c in  the  interior 
ofMn.  Then  there  exists  a unique  quadrature 

P rb 

u>i  fr(xi)  = fr(x)da(x),  r=l,...,/i,  (75) 

of  index  n /2,  where  Wi  > 0 and  x,-  € I,  i = 1, . . . , p.  In  particular,  if  n = 2m,  then  p — m and 

a < x\  < X2  < • • • < xm  < b\  (76) 

if  n — 2m  + 1 , then  p = m + \ and 

a < x i < X2  < ■ ■ • < xm+\  = b.  (77) 


Proof  The  Chebyshev  property  implies  that  there  exists  8 with  a < 8 < b such  that  /,  is  nonzero  on 
(a,8],i  = 1, . . . , n.  We  define  the  function  u on  / by  the  formula 


u(x)  = 


max  | /,(*)| 

/=1 n 

u(8) 


x € (a,  5], 
x € (8,  b], 
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and  observe  that  u is  continuous  and  positive  on  I and  integrable  on  [a,  b ] with  respect  to  a.  Now  we 
define  functions  g\, . . . , gn  on  [a,  b]  by  the  formula 


8i(x) 


M* ) 

u(x ) 

lint  gi(x) 


x e (a,  b], 
x — a, 


i = 1,...  ,n. 


The  system  g\ , . . . , gn  is  a Chebyshev  system  on  [< a , b]  and  is  integrable  with  respect  to  the  measure 
f*  u{t)da(t).  Theorem  4.2  therefore  is  applicable  and  assures  the  existence  of  exactly  two  quadra- 
tures of  index  nil  for  the  interval  [a,  b],  one  of  which  includes  the  point  x\  — a.  Our  assumption 
Xi  € (a,  b]  excludes  this  case  and  we  are  left  with  the  single  quadrature  presented  in  (76)  and  (77).  □ 


The  next  theorem,  which  is  the  principal  analytical  result  of  this  section,  follows  directly  from  The- 
orems 4.5  and  4.6.  The  existence  and  uniqueness  of  the  quadratures  defined  in  Section  3 follow  from 
it.  It  also  hints  at  the  existence  of  somewhat  more  general  quadratures,  for  singularities  of  the  form 
xy  log*  x,  but  we  do  not  evaluate  these  here. 

Theorem  4.7.  Suppose  yi, ...  ,Yj  are  distinct  real  numbers,  each  greater  than  —1,  and  n\,  . . . , n} 
are  positive  integers  with  ^ n,  = n.  For  sufficiently  large  a,  the  system  of  equations 

\n/r\ 

Y2  wmxmY'\og  xm  = (-l)*+1£w(- Yi,a),  k = 0, ...  , n,  - 1,  i = 1,...  ,j,  (78) 

m- 1 

has  a unique  solution  w i, ...  ,wn,X\,...  ,xn  satisfying  Wj  > 0,  i = 1, . . . , \n/T\,  and  0 < x\  < 
. . . < X[n/2-\  < a,  with  x [n/2]  = a if  n is  odd. 


5.  Computation  of  the  Nodes  and  Weights 

The  nodes  and  weights  of  the  quadratures  defined  in  Section  3 are  computed  by  numerically  solving 
the  nonlinear  systems  (22),  (24),  (36),  and  (39).  Conventional  techniques  for  this  problem  are  either 
overly  cumbersome  or  converge  too  slowly  to  be  practical.  Recently,  Ma,  Rokhlin,  and  Wandzura  [6] 
addressed  this  need  by  developing  a practical  numerical  algorithm  that  is  effective  in  a fairly  general 
setting.  They  construct  a simplified  Newton’s  method  and  combine  it  with  a continuation  (homotopy) 
method.  We  present  their  method  in  an  abbreviated  form  in  Section  5.2;  the  reader  is  referred  to  [6] 
for  more  detail. 

The  systems  for  regular  integrands,  however,  can  be  solved  even  more  simply,  as  we  see  next. 

5.1.  Regular  Integrands.  The  classical  theory  of  Gaussian  quadratures  for  polynomials,  summa- 
rized in  Section  2.4,  can  be  exploited  to  solve  (22)  and  (24).  In  particular,  suppose  that  po,  ■ ■ ■ , Pj 
are  the  orthogonal  polynomials  on  [0,  a],  given  by  the  recurrence  (13)  and  (14),  under  the  assumption 
that 

fa  r , w Br+i(a)  n „ . , 

/ x co(x)dx  = — — , r = 0, ...  ,2]  — 1. 

Jo  r + 1 

Then  the  roots  x\ , . . . , Xj  of  pj  and  corresponding  Christoffel  numbers  uq , . . . , Wj  satisfy  (22).  The 
polynomials  p\ , . . . , pj  can  be  calculated  in  symbolic  form;  their  coefficients  are  rational  if  a is  ratio- 
nal. The  roots  of  pj  can  be  computed  by  Newton  iteration  and  the  Christoffel  numbers  can  be  obtained 
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using  (16).  Similar  treatment  can  be  applied  to  the  system  (24)  containing  an  odd  number  of  equations, 
using  the  interval  [0,  a — 1],  under  the  assumption 

r Br+i(a) 

/ x x {x)dx  — , r = 0, . . . , 2j  — 2. 

Jo  r + 1 

The  Gauss-Radau  quadrature  is  computed  using  the  formula  ( 1 8)  for  the  modified  Christoffel  numbers. 
For  these  tasks  it  is  convenient  to  use  a software  system  that  can  manipulate  polynomials  with  full- 
precision  rational  coefficients.  The  author  implemented  code  for  these  computations  in  Pari/GP  [ 14]. 

It  should  be  noted  that  the  proposed  procedure  is  suitable  for  relatively  small  values  of  j (less  than, 
say,  20).  It  is  neither  very  efficient  nor  very  stable,  but  was  quite  adequate  for  our  purposes.  (Unlike 
the  situation  for  standard  Gaussian  quadratures,  where  the  number  of  nodes  depends  on  the  size  of  the 
problem,  here  j depends  only  on  the  desired  order  of  convergence.)  If  it  is  required  to  compute  the 
nodes  and  weights  of  (22)  or  (24)  for  large  j,  the  reader  may  consider  numerical  schemes  for  Gaussian 
quadrature  proposed  by  other  authors,  for  example,  that  of  Gautschi  [15]  or  Golub  and  Welsch  [16], 


5.2.  Singular  Integrands.  The  systems  (36)  and  (39)  for  singular  integrands  cannot  be  solved  us- 
ing methods  for  standard  Gaussian  quadrature,  since  the  nodes  to  be  computed  do  not  coincide  with 
the  roots  of  any  closely-related  orthogonal  polynomials.  We  employ  instead  the  algorithm  for  such 
systems  developed  by  Ma,  Rokhlin,  and  Wandzura  [6],  which  we  now  describe. 

A collection  of  2 n real-valued  continuously  differentiable  functions  f\ , . . . , f2n  defined  on  an  in- 
terval I — \a,  b]  is  an  Hermite  system  if 


fi{x\)  •• 

• f2n(*l)\ 

//(*  i) 

f2(x  1)  •• 

■ fin  (•*! ) 

f\  (*2) 

f2(x2)  ■ 

• fln(x2) 

det 

/1U2) 

f2(x2)  ■ 

■ flnfo) 

f\(xn) 

fliXn)  ■ 

fin  ( Xn  ) 

\f{(Xn) 

f2(xn)  • 

• flniXn)/ 

for  any  choice  of  distinct  x\, . . . 

, xn  on  /. 

An  Hermite  system  that 

7^0 


(79) 


extended  Hermite  system.  The  following  theorem  is  a direct  consequence  of  the  definition;  the  proofs 
of  the  subsequent  two  theorems  are  contained  in  [6]. 


Theorem  5.1.  Suppose  that  the  functions  f\ , . . . , f2n  constitute  an  Hermite  system  on  the  interval 
[a,  b ] and  jcj,  . . . , xn  are  n distinct  points  on  [a,  b].  Then  there  exist  unique  coefficients  otjj , , i = 

1, ...,«,  7 = 1, ...  , 2n,  such  that 


Oi(xk)  = 0 

rh(xk)  = 8ik 

(80) 

(y'(xk)  = 8ik 

h'i(xk)  = 0 

(81) 

"t 

II 

II 

. , n,  where  the  functions  Oj,  r/j  are  defined  by  the  formulae 

In 

°i(x)  = Yla,j 

7 = 1 

In 

(82) 

rii(x)  = X]  f‘J  fjW 

7=1 

(83) 

for  i = 1 , . . . , n. 
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Theorem  5.2  ([6]  p.  979).  Suppose  that  the  functions  fi, ...  , fin  are  an  Hermite  system  on  [a,  b ]. 
Suppose  further  that  S C [a,  b]n  is  the  set  of  points  with  distinct  coordinates  X\, . . . ,xn.  Suppose 
finally  that  the  mapping  F : S -»  K."  is  defined  by  the  formula 


F(x\, ...  ,xn)=  yj  oi(x)a>(x)dx, ...  , J an(x)  <o(x)  dx^j . 


(84) 


with  the  functions  a\, . . . , an  defined  by  (80)  through  (82).  Then  x\ , . . . , xn  are  the  Gaussian  nodes 
for  the  system  of  functions  f\ , . . . , fin  with  respect  to  the  weight  co  if  and  only  if  F(x\, ...  , xn)  — 0. 


Theorem  5.3  ([6]  p.  983).  Suppose  that  the  functions  fi, . . . , fin  are  an  extended  Hermite  system 
on  [a,  b],  S C [a,  b]n  is  the  set  of  points  with  distinct  coordinates,  and  the  mapping  G : S — ► M"  is 
defined  by  the  formula 


G(x i, ...  ,xn)  = 


fa  ai(x)co(x)dx 

fa  m(x)a)(x)dx 


. . . , xn  + 


fa  crn(x)  (o(x)  dx 
fa  hn(x)co(x)dx 


(85) 


with  the  functions  cr\ , . . . ,on  and  rj  i , . . . , tjn  defined  by  (80)  through  (83).  Suppose  further  that  fi  € 
C3((a,  b))  for  i = l, ...  ,2  n and  the  function  F is  defined  by  (84).  Suppose  finally  that  x*  is  the 
unique  zero  of  F,  that  Xo  is  an  arbitrary  point  of  S,  and  that  the  sequence  Xj , Xj, . . . is  defined  by  the 
formula 


x,-+i  =G(Xi),  i = 0,1 


(86) 


Then  there  exists  € > 0 and  a > 0 such  that  the  sequence  Xi,  Xi, . . . generated  by  (86)  converges  to 
x*  and 

l|x;+i  - x*||  < or  ||x,- - x*||2  (87) 


for  any  initial  point  xq  such  that  ||xq  — x*||  < e. 


The  key  feature  of  this  theorem  is  the  quadratic  convergence  indicated  by  (87).  The  solution  x* 
is  obtained  by  an  iterative  procedure;  each  step  consists  of  computing  the  coefficients  that  determine 
cti, ...  ,csn  and  rji, . . . , qn  by  inverting  the  matrix  of  (79),  then  computing  the  integrals  that  define 
G by  taking  linear  combinations,  using  these  coefficients,  of  integrals  of  fi, . . . , f Theorem  5.3 
ensures  that  with  appropriate  choice  of  starting  value  xo,  convergence  is  rapid  and  certain. 

For  quadrature  nodes  x*  = (*i, . . . , xn ),  the  quadrature  weights  are  given  by  the  integrals  of  77,, 
namely 


(ici, . . . ,wn)  = 


qi(x)  a>(x)  dx. 


T}\(x)a)(x)  dx  ) . 


(88) 


We  note  that  Theorems  5.1  through  5.3  concern  Gaussian  quadratures  with  n nodes  and  weights  to 
integrate  2 n functions  on  the  interval  [a , b ] exactly.  For  Gauss-Radau  quadratures,  in  which  node  xn  = 
b is  fixed  and  2n  — 1 functions  are  integrated  exactly,  only  a slight  change  is  required.  In  particular, 
functions  o\ , . . . , cr„_t  (without  an)  and  r/\, . . . , rjn  are  defined  as  before  by  (80)  through  (83),  except 
that  the  summations  in  (82)  and  (83)  exclude  fin.  Their  coefficients  aq,  i = 1, . . . , n — 1,  j = 
1, . . . ,2  n — 1,  and  fiq,  i = 1 , ...  ,n,  j = 1, . . . , 2n  — 1,  are  obtained  by  inverting  the  matrix 
which  results  from  removing  the  last  row  and  column  from  the  matrix  of  (79).  The  revised  mapping 
G : S — > M”-1  has  n — 1 components  defined  as  the  first  n — 1 components  in  (85).  Finally,  as  before, 
the  quadrature  weights  are  given  by  (88). 


HYBRID  GAUSS-TRAPEZOIDAL  QUADRATURE  RULES 


21 


In  order  to  obtain  a sufficiently  good  starting  estimate  for  the  solution  of  F(x)  = 0,  a continuation 
procedure  can  be  used,  as  outlined  in  the  following  theorem. 

Theorem  5.4  (see,  for  example,  [6]  p.  975).  Suppose  that  F : [0,  1]  x [a,  b]n  ->  R"  is  a function 
with  a unique  solution  x,  to  the  equation  F(f , x)  = 0 for  all  t € [0,  1],  suppose  that  x,  is  a continuous 
function  oft,  and  suppose  that  F(1 , x)  = Fix).  Finally,  suppose  Xo  is  given  and  for  some  8 > 0 there 
is  a procedure  P to  compute  xt  for  t G [0,  1],  given  an  estimate  x,  with  |x,  — x,|  < 8.  Then  there 
exists  a positive  integer  m such  that  the  following  procedure  can  be  used  to  compute  the  solution  of 
F(x)  = 0: 


For  i = 1, . . . , m,  use  P to  compute  x,/m,  given  the  estimate  X(,_  \)/m. 

The  required  solution  of  Fix)  = 0 isx\. 

More  typically,  of  course,  8 and  any  bound  on  |x,+e  — x,  \/e  depend  on  t and  in  a practical  implemen- 
tation the  step  size  is  chosen  adaptively. 

To  compute  the  solutions  of  (36)  and  (39),  it  is  effective  to  use  a continuation  procedure  with  respect 
to  both  j and  a.  Solutions  for  the  first  few  values  of  j are  readily  obtained  without  requiring  good 
initial  estimates.  Given  a solution  of  (36)  for  the  interval  [0,  a ] with  nodes  ui, ...  ,uj  and  weights 
vi, . . . , Vj,  we  choose  an  initial  estimate  fij, . . . , uJ+\,  v i, . . . , vj+ \ for  j + 1 and  the  interval  [0,  a- f 
1]  defined  by  the  formulae 

(Ui  i = l,...,  j, 
a i = j + 1, 

This  choice  exactly  satisfies  the  equations 

r.  Hj  Viy+r  = -C(-y  - r,  a + 1), 

;=i 

- ~ r Br+ 1 ia  + 1) 

Ui  Vi  = — , 

r + 1 

as  follows  immediately  from  the  difference  formula  (4)  for  Bn  and  (11)  for  f,  but  the  corresponding 
equations  for  r = j are  not  satisfied.  Those  equations  are  approximately  satisfied,  however,  and  we 
can  start  with  the  actual  values  of  the  sums  for  r = j as  the  required  values.  These  are  then  var- 
ied continuously,  obtaining  the  corresponding  solutions,  until  they  coincide  with  the  intended  values 
— f (— y — j,a- f 1)  and  Bj+\ia  + 1 )/ij  + 1).  This  procedure  can  be  used  without  alteration  for  (39). 

Once  the  solution  for  j + 1 and  the  interval  [0,  a + 1]  is  obtained,  a can  be  continuously  varied,  in 
a continuation  procedure,  to  obtain  solutions  for  different  intervals. 

6.  Numerical  Examples 

The  procedures  described  in  Section  5 were  implemented  in  Pari/GP  [14]  for  both  the  regular  cases 
and  the  singular  cases.  The  matrix  in  (79),  which  must  be  inverted,  is  very  poorly  conditioned  for  many 
choices  of  n,  x\, . . . , xn,  and  f\, . . . , fzn.  This  difficulty  was  met  by  using  the  extended  precision 
capability  of  Pari/GP. 


L 


r = 0,  1, . . . , j — 1, 
r = 0,  1, — 1, 


(90) 


Vi  = 


Vj  ( = 1,...,  j, 

1 1=7  + 1. 


(89) 
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Table  1.  The  minimum  value  of  a,  considered  as  a function  of  j,  such  that  the 
point  [B\{a)/\, , Bj(a)/j ) is  in  the  moment  space  !M2]  of  the  polynomials 
1,  x, . . . , x2j~l  on  the  interval  [0,  a].  The  moment  space  is  defined  at  (52). 


j 

min  a 

j 

min  a 

j 

min  a 

j 

min  a 

1 

0.78868 

5 

3.96696 

9 

7.21081 

13 

10.47885 

2 

1.57085 

6 

4.77448 

10 

8.02618 

14 

11.29815 

3 

2.36347 

7 

5.58463 

11 

8.84274 

15 

12.11815 

4 

3.16288 

8 

6.39687 

12 

9.66035 

16 

12.93878 

6.1.  Nodes  and  Weights.  The  nodes  and  weights  of  (22),  (24),  (36),  and  (39)  that  determine  the 
quadratures  of  Section  3 were  computed  for  a range  of  values  of  the  parameter  j.  For  each  choice  of 
j,  a was  chosen,  by  experiment,  to  be  the  smallest  integer  leading  to  positive  nodes  and  weights  (see 
Theorem  4.7).  For  the  regular  case  (22),  the  characterization  in  Theorem  4.4  was  used  to  determine  the 

minimum  value  of  a € R for  j — 1, . . . ,16,  that  satisfies  c = [B\(a)/\ Bj(a)/j ) e M2  j-  In 

particular,  we  obtained  the  minimum  value  of  a such  that  the  quadratic  forms  in  (62)  are  nonnegative 
definite.  This  determination  was  made  by  calculating  the  determinant  of  each  corresponding  matrix 
symbolically  and  solving  for  the  largest  root  of  the  resulting  polynomial  in  a.  These  values  are  given 
in  Table  1.  This  evidence  suggests  that  lirn^oc  j~l  min  a exists  and  is  roughly  5/6,  meaning  that  the 
number  of  trapezoidal  nodes  displaced  is  less  than  the  number  of  Gaussian  nodes  replacing  them,  for 
quadrature  rules  of  all  orders.  This  relationship  also  appears  to  hold,  to  an  even  greater  extent,  for  the 
singular  cases. 

The  values  of  selected  nodes  and  weights,  for  the  regular  cases  and  for  singularities  1/2  and  log*, 

are  presented  in  an  appendix.  Of  particular  simplicity  are  the  first  two  rules  for  regular  integrands. 


T'A(f)=h 


+ /(l  — h)  + —f  (1  — h/6) 


(91) 


where  h = 1 /(n  + 1 ),  for  n =0,  1 , . . . , and 


T~2(f)=h 


!/<A/5)  + S/(A)  + /<2/,)  + 


+/( i - 2*>  + - *)  + § /O  - A/5) 


(92) 


where  h = 1 /(n  + 3),  for  n =0,  1 These  rules  are  of  third  and  fourth  order  convergence,  respec- 

tively. The  first  is  noteworthy  for  having  the  same  weights  as,  but  higher  order  than,  the  trapezoidal 
rule;  the  second  has  asymptotic  error  1 /4  that  of  Simpson’s  rule  with  the  same  number  of  nodes. 

The  lowest-order  rule  presented  for  logarithmic  singularities. 


’l.i.u 


(g)  = h 


-g(h/(2it))  + g(h)  + 


+ g(l-h)  + -g(l) 


(93) 


approximates  fQl  g(x)  dx  with  error  of  order  0(h 2 log h)  for  g(x)  = (p(x ) logx  + \ Js(x),  provided  0 
and  0-  are  regular  functions  on  [0,  1],  The  corresponding  rule  for  the  singularity  x_1/2  is  not  quite  as 
simple,  for  there  the  quantity  2jt  in  (93)  is  replaced  with  4£(l/2)2. 
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6.2.  Quadrature  Performance.  To  demonstrate  the  performance  of  the  quadrature  rules,  they  were 
used  in  a Fortran  routine  (with  real*  8 arithmetic)  to  numerically  compute  the  integrals 

l 

[cos(200x)  5(r)  + cos(200r  -f  .3)]  dx  (94) 

for  the  functions  s(r)  = 0,  s(x ) = r_1/2,  and  .y(jt)  = log*.  These  integrals  were  also  obtained 
analytically  and  the  relative  error  of  the  quadratures  was  computed.  The  numerical  integrations  were 
computed  for  various  orders  of  quadrature  and  various  numbers  of  nodes.  Minimum  sampling  was 
taken  to  be  two  points  per  period  of  the  cosine  (i.e.,  200/7T  ~ 63.7  quadrature  nodes).  The  accuracies 
are  then  compared  for  various  degrees  of  oversampling.  The  quadrature  errors  are  listed  in  Tables  2 
through  4 and  plotted,  as  a function  of  oversampling  factor,  in  Figure  1 . We  note  that  the  graphs 


Table  2.  Relative  errors  in  the  computation  of  the  integral  in  (94),  for  the  regular 
case  s(x)  = 0.  Quadrature  rules  with  convergence  of  order  2,  4,  8,  16,  and  32  were 
used  with  various  numbers  m = n + 2j  of  nodes.  Here  / = m/( 200/tt)  is  the 
oversampling  factor. 


m 

/ 

2 

4 

8 

16 

32 

70 

1.10 

0.622D+00 

0.114D— 01 

0.382D— 02 

0.170D— 05 

0.234D— 10 

80 

1.26 

0.488D+00 

0.938D— 02 

0.184D— 02 

0.354D— 06 

0.115D— 11 

90 

1.41 

0.391D+00 

0.744D— 02 

0.934D— 03 

0.841D— 07 

0.720D— 13 

100 

1.57 

0.321D+00 

0.584D— 02 

0.498D— 03 

0.223D— 07 

0.192D— 14 

115 

1.81 

0.246D+00 

0.408D— 02 

0.21  ID— 03 

0.365D— 08 

0.331D— 14 

130 

2.04 

0.194D+00 

0.289D— 02 

0.964D— 04 

0.715D— 09 

0.331D— 14 

145 

2.28 

0.157D+00 

0.209D— 02 

0.472D— 04 

0.162D— 09 

0.471D— 14 

160 

2.51 

0.129D+00 

0.154D— 02 

0.245D— 04 

0.415D— 10 

0.262D— 14 

180 

2.83 

0.102D+00 

0.106D— 02 

0. 110D— 04 

0.794D— 11 

0.471D— 14 

200 

3.14 

0.832D— 01 

0.747D— 03 

0.531D— 05 

0.177D— 11 

0.33  ID- 14 

230 

3.61 

0.631D— 01 

0.465D— 03 

0.199D— 05 

0.235D— 12 

0.523D— 15 

260 

4.08 

0.495D— 01 

0.303D— 03 

0.828D— 06 

0.375D— 13 

0.384D— 14 

Table  3.  Relative  errors  for  the  singular  case  s(x)  = x 1/2,  for  various  numbers 
m = n + j + k of  nodes  and  orders  of  convergence. 


m 

/ 

2 

4 

8 

16 

70 

1.10 

0.692D— 01 

0.519D— 01 

0.850D— 02 

0.163D— 03 

80 

1.26 

0.925D— 01 

0.258D— 01 

0.260D— 02 

0.578D— 05 

90 

1.41 

0.921D— 01 

0.133D— 01 

0.698D— 03 

0.667D— 06 

100 

1.57 

0.838D— 01 

0.717D— 02 

0.146D— 03 

0.277D— 06 

115 

1.81 

0.686D— 01 

0.307D— 02 

0.201D— 04 

0.360D— 07 

130 

2.04 

0.550D— 01 

0.144D— 02 

0.269D— 04 

0.437D— 08 

145 

2.28 

0.441D— 01 

0.724D— 03 

0.171D— 04 

0.557D— 09 

160 

2.51 

0.357D— 01 

0.389D— 03 

0.964D— 05 

0.733D— 10 

180 

2.83 

0.273D— 01 

0.186D— 03 

0.440D— 05 

0.408D— 11 

200 

3.14 

0.212D— 01 

0.976D— 04 

0.207D— 05 

0.218D— 12 

230 

3.61 

0.151D— 01 

0.427D— 04 

0.724D— 06 

0.130D— 12 

260 

4.08 

0. 110D— 01 

0.215D— 04 

0.280D— 06 

0.201D— 13 
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Table  4.  Relative  errors  for  the  singular  case  s(x)  = log*.  Here  the  error  is  of 
order  0(hl  log  h),  where  / is  shown. 


m 

/ 

2 

4 

8 

16 

70 

1.10 

0.369D+00 

0.217D— 01 

0.354D— 01 

0.243D— 03 

80 

1.26 

0.271D+00 

0.238D— 02 

0.328D— 02 

0.487D— 04 

90 

1.41 

0.206D+00 

0.765D— 02 

0.707D— 03 

0.394D— 05 

100 

1.57 

0.162D+00 

0.768D— 02 

0.687D— 03 

0.121D— 05 

115 

1.81 

0.117D+00 

0.576D— 02 

0.291D— 03 

0.886D— 07 

130 

2.04 

0.882D— 01 

0.398D— 02 

0.120D— 03 

0.903D— 08 

145 

2.28 

0.687D— 01 

0.272D— 02 

0.548D— 04 

0.123D— 08 

160 

2.51 

0.549D— 01 

0.188D-02 

0.272D— 04 

0.177D— 09 

180 

2.83 

0.421D— 01 

0.119D— 02 

0.118D-04 

0.965D— 11 

200 

3.14 

0.332D— 01 

0.774D— 03 

0.550D— 05 

0.956D— 12 

230 

3.61 

0.243D— 01 

0.433D— 03 

0.196D— 05 

0.398D— 12 

260 

4.08 

0.185D— 01 

0.258D— 03 

0.778D— 06 

0.106D— 12 

s(x)=0 

s(x)=xA(- 

1/2) 

s(x)=log(x 

1 2 3 4 

Oversampling  (f) 


Oversampling  (f) 


Oversampling  (f) 


FIGURE  1 . The  relative  errors  of  the  quadratures,  shown  in  Tables  2 through  4,  are 
plotted  using  logarithmic  scaling  on  both  axes. 


are  nearly  straight  lines  (until  the  limit  of  machine  precision  is  reached),  as  predicted  from  the  theo- 
retical convergence  rates.  We  remark  also  that  excellent  accuracy  is  attained  for  even  quite  modest 
oversampling  when  quadratures  with  high-order  convergence  are  employed.  For  problems  where  the 
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Table  5.  Relative  errors  in  the  computation  of  the  integral  in  (95),  for  £ = 1.  Quad- 
rature rules  defined  in  Corollary  3.10  with  j = 1,  2,  4,  8,  and  16  were  used  with 
various  numbers  m of  nodes. 


m 

1 

2 

4 

8 

16 

70 

0.999D+00 

0.400D+00 

0.305D+00 

0.180D+00 

0.104D— 01 

80 

0.304D+00 

0.200D— 01 

0.247D— 01 

0.238D— 02 

0.474D— 03 

90 

0.113D+00 

0.217D— 01 

0.136D— 02 

0.383D— 03 

0.866D— 05 

100 

0.273D— 01 

0.137D— 01 

0.137D— 02 

0.440D— 04 

0.900D— 06 

115 

0.210D— 01 

0.247D— 02 

0.137D— 03 

0.573D— 05 

0.331D— 07 

130 

0.228D— 01 

0.107D— 02 

0.632D— 04 

0.423D— 06 

0.107D— 09 

145 

0.118D— 01 

0.115D— 02 

0.305D— 04 

0.196D— 06 

0.490D— 10 

160 

0.21 2D— 02 

0.521D— 03 

0.307D— 05 

0.430D— 07 

0.216D— 09 

180 

0.625D— 02 

0.824D— 04 

0.451D— 05 

0.101D— 07 

0.166D— 11 

200 

0.558D— 02 

0.206D— 03 

0.184D— 05 

0.291D— 08 

0.867D— 12 

230 

0.863D— 03 

0.676D— 04 

0.463D— 06 

0.626D— 09 

0.586D— 13 

260 

0.266D— 02 

0.433D— 04 

0.291D— 06 

0.596D— 10 

0.346D— 14 

number  of  quadrature  nodes  is  the  major  cost  factor,  therefore,  one  may  benefit  by  using  the  high-order 
quadratures  even  for  modest  accuracy  requirements. 

We  test  the  quadratures  for  improper  integrals  by  numerically  computing  for  £ = 1 the  integral 

/ e~ix*  J"  : dx  — —2ni  //(£)  Y"  (r  + l)eir^,  (95) 

•'-0°  10-*+r+Z  r=— 10 

where  H is  the  Heaviside  step  function.  The  integrand  is  oscillatory  and  decays  like  a:-1  as  x ±oo. 
The  quadratures  defined  in  Corollary  3.10  are  employed,  for  which  the  integral  is  split  into  a regular 
integral  on  a finite  interval,  chosen  here  to  be  [-5-^/4,  5y/m/4],  where  m is  the  total  number  of 
quadrature  nodes,  and  two  improper  integrals  in  the  imaginary  direction,  using  Laguerre  quadratures. 
The  quadrature  errors  are  shown  in  Table  5. 

7.  Applications  and  Summary 

The  chief  motivation  for  the  hybrid  Gauss-trapezoidal  quadrature  rules  is  the  accurate  computation 
of  integral  operators.  We  define  an  integral  operator  A by  the  formula 

(Af)(x)  = J K(x,  y)  f (y)dy, 

where  T is  a regular,  simple  closed  curve  in  the  complex  plane,  the  function  / is  regular  on  T,  and 
the  kernel  K : C x C ->  C is  a regular  function  of  its  arguments,  except  where  they  coincide;  we 
assume 

K(x,  y ) = <t>(x,  y)s(|*  - y\)  + x/s(x,  y),  (96) 

with  (f)  and  i/r  regular  on  T x T and  5 regular  on  (0,  oo),  with  an  integrable  singularity  at  0.  A large 
variety  of  problems  of  classical  physics  can  be  formulated  as  integral  equations  that  involve  such  op- 
erators. When  the  operator  occurs  in  an  integral  equation 

f{x)  + ( Af){x ) = g(x). 


x € r, 


(97) 
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some  choice  of  discretization  must  be  used  to  reduce  the  problem  to  a finite-dimensional  one  for  nu- 
merical solution.  In  the  Nystrom  method  the  integrals  are  replaced  by  quadratures  to  yield  the  finite 
system  of  equations 

m 

f(xi)  + ^2WU  f(Xj)  = g(Xi),  i = l, . . . , m.  (98) 

7 = 1 

This  linear  system  can  be  solved  for  f(x\),...,f(xm)bya.  variety  of  techniques.  The  particular 
choice  of  and  wtj  for  i,  j = l, ...  ,m  determines  the  order  of  convergence  (and  therefore  effi- 
ciency) of  the  method. 

For  a curve  parametrization  v : [0,  1]  — > T,  such  as  scaled  arc  length,  the  operator  A becomes 

(Af)(v(t))=  [ K(v(t),  v(r))  /(v(r))  v'(r)  dr. 

Jo 

It  is  convenient  to  use  a uniform  discretization  1/m,  2/m, . , . , 1 in  t and  r,  so  jc,  = v(i / m ),  i = 
1, . . . , m.  How  then  is  ic,7  determined?  We  assume  for  the  moment  that  / is  available  at  locations 
other  than  x \ , . . . ,xm.  Continuing  v periodically  with  period  1 , and  using  the  Gauss-trapezoidal  quad- 
ratures, we  obtain 


(Af)(v(i/m)) 


-r 

J i/m 


1 +i/m 


K(v(i/m),  v(r))  f(v( r)  v'(z)dz 


1 J j n-l  j J 

~ UkCTj/mivk/m)  + — <Jiim(a/m  +k/m)  + — V uk  ai/m(l  - vk/m)  (99) 
m mtJo 

for  / = 1 , . . . , m,  where  oa  : [0,  1]  — > C is  defined  by  the  formula 

cra(r)  = K(v(a),  v{ot  + r))  f(a  + r)  v'(ct  + r)  (100) 


and  m = n + 2a  — \ and  u \, . . . ,Uj,v i , . . . , vj  are  determined  for  the  singularity  s of  K.  Provided 
that  the  periodic  continuation  of  v is  sufficiently  regular,  the  quadrature  will  converge  to  the  integral 
with  order  greater  than  j as  m — > oo,  for  / = 1 , . . . , m.  We  relax  the  restriction  that  / be  available 
outside  x\ , . . . , xm  by  using  local  Lagrange  interpolation  of  order  j + 1 for  equispaced  nodes. 


f(v (t))  % + r/m ))  Irimx  - i ), 


r= 0 


where  i = [mz  — ( j — 1)/2J  and 


lr(x)=  FT  , r = 0, . . . , j. 

sAs*rr~S 


(101) 


Now  Wij  for  i,  j = l, ...  ,m  is  determined  by  combining  (97)  through  (101).  The  computation  of 
all  m 2 coefficients  requires  m(m  +2  j —2a  + 1)  evaluations  of  the  kernel  K and  therefore  order 
0(m2)  operations.  This  cost  can  often  be  substantially  reduced  using  techniques  that  exploit  kernel 
smoothness  (see,  for  example,  [17],  [18]). 

A slightly  different  application  of  the  quadratures  is  the  computation  of  Fourier  transforms  of  func- 
tions that  fail  to  satisfy  the  assumptions  usually  made  when  using  the  discrete  Fourier  transform.  In 
particular,  if  a function  decays  slowly  for  large  argument  or  is  compactly  supported  and  singular  at 
the  ends  of  the  support  interval,  these  quadratures  can  be  used  to  compute  its  Fourier  transform.  One 


HYBRID  GAUSS-TRAPEZOIDAL  QUADRATURE  RULES 


27 


example  of  such  a function  is  that  in  (95).  Since  most  of  the  nodes  in  these  quadratures  are  equispaced, 
with  function  values  given  equal  weight,  the  fast  Fourier  transform  can  be  used  to  do  the  bulk  of  the 
computations;  the  overall  complexity  is  0(n  logrt),  where  n is  the  number  of  Fourier  coefficients  to 
be  computed. 

Other  applications  may  include  the  representation  of  functions  for  solving  ordinary  or  partial  differ- 
ential equations,  when  high-order  methods  are  required.  In  addition,  an  extension  of  these  quadratures 
to  integrals  on  surfaces  is  under  study. 

To  summarize,  the  characteristics  of  the  hybrid  Gauss-trapezoidal  quadrature  rules  include 

• Arbitrary  order  convergence  for  regular  functions  or  functions  with  known  singularities  of  power 
or  logarithmic  type, 

• Positive  weights, 

• Most  nodes  equispaced  and  most  weights  constant,  and 

• Invariant  nodes  and  weights  (aside  from  scaling)  as  the  problem  size  increases. 

The  primary  disadvantage  of  the  quadrature  rules,  shared  with  other  Gaussian  quadratures  but  exacer- 
bated here  by  poor  conditioning,  is  that  the  computation  of  the  nodes  and  weights  is  not  trivial.  Nev- 
ertheless, tabulation  of  nodes  and  weights  for  a given  order  of  convergence  allows  this  issue  to  be 
avoided  in  the  construction  of  high-order,  general  purpose  quadrature  routines. 
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Appendix.  Tables  of  Quadrature  Nodes  and  Weights 

Table  6.  The  nodes  and  weights  for  the  quadrature  rule  Tila{f)  = 

h 5Z/=i  w'  f (xih)  + h o flflh  + ih)  + h EL  w‘  /O  - *«*).  with 
h — (n  + 2a  — l)-1,  for  several  choices  of  j and  corresponding  minimum 
integer  a.  For  / a regular  function,  T„a{f)  converges  to  /0'  f(x)dx  as  n — ► oo 
with  convergence  of  order  O. 


o 

a 

Xi 

Wi 

3 

1 

1.66666  66666  66667D-01 

5.00000  00000  00000D— 01 

4 

2 

2.00000  00000  00000D— 01 

1 .00000  00000  00000D+00 

5.20833  33333  33333D-01 
9.79166  66666  66667D-01 

5 

2 

2.24578  49798  12614D-01 
1.01371  93743  59164D+00 

5.54078  16436  06372D-01 
9.45921  83563  93628D-01 

6 

3 

2.25099  10426  10971D-01 
1.01426  90609  87992D+00 
2.00000  00000  OOOOOD+OO 

5.54972  43271  64180D-01 
9.45131  74118  45473D-01 
9.99895  82609  90347D-01 

7 

3 

2.18054  06725  43505D-01 
1.00118  18730  31216D+00 
1.99758  05264  18033D+00 

5.40808  89672  08193D-01 
9.51661  50458  23566D-01 
1.00752  95986  96824D+00 

8 

4 

2.08764  74220  32129D-01 
9.78608  73737  14483D-01 
1.98954  13865  79751D+00 
3.00000  00000  OOOOOD-fOO 

5.20798  82772  46498D-01 
9.53503  80185  55888D-01 
1.02487  16264  02471D+00 
1.00082  57440  17291D+00 

12 

5 

7.02395  54616  21939D-02 
4.31229  78572  27970D-01 
1.11775  27345  18115D+00 
2.01734  37245  72518D+00 
3.00083  78428  47590D+00 
4.00000  00000  00000D+00 

1.92231  59778  43698D-01 
5.34839  95305  14687D-01 
8.17020  94424  88760D-01 
9.59211  15214  45966D-01 
9.96714  34080  44999D-01 
9.99982  01196  61890D-01 

16 

7 

9.91933  78414  51028D-02 
5.07659  26696  45529D-01 
1.18497  29258  27278D+00 
2.04749  34671  34072D+00 
3.00716  89118  693 10D+00 
4.00047  49967  76184D+00 

5.00000  78790  22339D+00 

6.00000  00000  00000D+00 

2.52819  89287  66921D-01 
5.55015  82301  59486D-01 
7.85232  14536  15224D-01 
9.24591  56738  76714D-01 
9.83935  02004  45296D-01 
9.98446  34484  13151D-01 
9.99959  23784  64547D-01 
9.99999  96862  58662D-01 

20 

9 

9.20920  04462  33291D-02 
4.75202  19477  58861D-01 
1.12468  79458  44539D+00 

1 .97738  73856  42367D+00 
2.95384  89578  22108D+00 
3.97613  67860  48776D+00 
4.99435  42819  79877D+00 
5.99946  95393  35291D+00 
6.99998  67048  74333D+00 
8.00000  00000  00000D+00 

2.35183  61446  43984D-01 
5.24882  05090  85946D-01 
7.63402  64098  69887D-01 
9.28471  13366  58351D-01 
1.01096  98865  87741D+00 
1.02495  97253  11073D+00 
1.01051  75346  39652D+00 
1.00155  15957  97932D+00 
1.00006  16817  94 188D+00 
1.00000  01358  43597D+00 
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Table  6.  (continued) 


o 

24 


28 


32 


a x, 

10  6.0010647314  74805D-02 
3.14968  50162  29433D-01 
7.66450  82405  18316D-01 
1.39668  57813  42510D+00 
2.17519  59032  06602D-I-00 
3.06232  05758  80355D+00 
4.01644  09887  92476D+00 
5.00287  20642  75734D+00 
6.00028  54533  10164D+00 
7.00001  29649  62529D+00 

8.00000  01755  54469D+00 

9.00000  00000  00000D+00 
12  6.23436  05331  94102D-02 

3.25028  67217  02614D-01 
7.83735  07942  82182D-01 
1.41567  31126  16924D+00 
2.18989  42500  6131 3D+00 
3.07005  38774  83040D+00 
4.01861  37562  18047D+00 
5.00270  59020  35397D+00 
5.99992  97418  10400D+00 
6.99990  47208  46024D+00 

7.99998  68948  43540D+00 

8.99999  93733  80393D+00 

9.99999  99920  0291 1D+00 
1 . 1 0000  00000  00000D+0 1 

14  5.89955  06143  25259D-02 
3.08275  70622  27814D-01 
7.46370  72530  79130D-01 
1.35599  37264  94664D+00 
2.11294  32173  46336D+00 
2.98724  14965  45946D+00 
3.94479  89209  61 176D+00 
4.95026  92028  42798D+00 
5.97212  30431  17706D+00 
6.98978  35581  37742D+00 
7.99767  30195  12965D+00 
8.99969  49327  47039D+00 
9.99997  92252  11805D+00 
1.09999  99382  66130D+01 
1.19999  99994  62073D+01 
1 .30000  00000  OOOOOD+Ol 


t m 

1.53893  21045  18340D-01 
3.55105  81285  59424D-01 
5.44920  00362  80007D-01 
7.10407  84977  15549D-01 
8.39878  09402  53654D-01 
9.27276  79508  90611D-01 
9.75060  56973  71132D-01 
9.94262  96508  23470D-01 
9.99242  17784  21898D-01 
9.99953  43707  86161D-01 

9.99999  08549  12925D-01 

9.99999  99894  66828D-01 
1.59597  52797  34157D-01 
3.63704  60281  93864D-01 
5.49875  31772  97441D-01 
7.08798  67920  86956D-01 
8.33517  22755  01 195D-01 
9.20444  65106  08518D-01 
9.71088  17765  52090D-01 
9.93329  65785  55239D-01 
9.99475  90879  10050D-01 
1.00013  30302  54421D+00 
1.00003  29150  11460D+00 

1.00000  22616  53775D+00 

1 .00000  00423  93520D+00 

1 .00000  00000  42872D+00 
1.51107  60238  74 179D-01 
3.45939  59211  69090D-01 
5.27350  28051  46873D-01 
6.87844  40945  43021D-01 
8.21031  91400  34114D-01 
9.21838  28755  15803D-01 
9.87302  74875  53060D-01 
1.01825  19134  41 155D+00 
1.02193  34303  49293D+00 
1.01256  79834  13513D+00 
1.00405  22895  54521D+00 
1.00071  34133  44501D+00 
1.00006  36183  02950D+00 

1.00000  24863  85216D+00 

1 .00000  00304  04477D+00 

1 .00000  00000  20760D+00 
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Table  7.  The  nodes  v\, . . . , Vj  and  weights  u\, . . . , uj  for  the  quadrature  rule 
SJnkab(g)  = h E/=i  Hi  g{Vih)  + h YHZg  g(ah  + ih)  + h £j=1  w,  g(l  - xth),  with 
h = (n+a+b  — l)-1,  for  g(.x)  = x~l/2<p(x)+\J/(x),  with  </>  and  xjr  regular  functions. 
The  nodes  x\, . . . , Xk  and  weights  w\ , . . . , are  found  in  Table  6. 


o 

a 

Vi 

Ui 

1.5 

1 

1.17225  85713  93266D-01 

5.00000  00000  00000D— 01 

2.0 

2 

9.25211  27154  21378D-02 

1 .00000  00000  00000D— 00 

4.19807  96252  66162D-01 
1.08019  20374  73384D+00 

2.5 

2 

6.02387  37964  08450D-02 
8.7807040506  76215D-01 

2.85843  99904  20468D-01 
1.21415  60009  57953D+00 

3.0 

2 

7.26297  84134  70474D-03 
2.24632  55125  2 1893D-01 

1 .00000  00000  OOOOOD+OO 

3.90763  87675  31813D-02 
4.87348  40566  46474D-01 
9.73575  20666  00344D-01 

3.5 

2 

1.28236  89094  58828D-02 
2.69428  63467  92474D-01 
1.0184145237  86358D+00 

6.36399  66631  05925D-02 
5.07743  45780  43636D-01 
9.28616  57556  45772D-01 

4.0 

3 

1.18924  24340  21285D-02 
2.57822  04347  38662D-01 
1.00775  00645  85281D+00 
2.00000  00000  00000D+00 

5.92721  50356  16424D-02 
4.95598  17403  06228D-01 
9.42713  12906  28058D-01 
1.00241  65465  50407D+00 

6.0 

4 

3.31792  59426  99451D-03 
8.28301  97052  96352D-02 
4.13609  49257  26231D-01 
1.08874  43736  88402D+00 
2.00648  21018  52379D+00 
3.00000  00000  00000D+00 

1.68178  09298  83469D-02 
1.75524  44045  44475D-01 
5.03935  05038  58001D-01 
8.26624  13396  80867D-01 
9.77306  58489  81277D-01 
9.99791  98099  47032D-01 

8.0 

5 

1.21413  06065  23435D-03 
3.22395  27000  27058D-02 
1.79093  53836  49920D-01 
5.43766  38052  4463 ID-01 
1.17611  66283  96759D+00 
2.03184  82107  16014D+00 
3.00196  12256  908 12D+00 
4.00000  00000  00000D+00 

6.19984  48842  97793D-03 
7.10628  67917  20044D-02 
2.40893  01044  10471D-01 
4.97592  92636  68960D-01 
7.59244  65404  41226D-01 
9.32244  63996  14420D-01 
9.92817  14381  60095D-01 
9.99944  91256  89846D-01 

10.0 

6 

1.74586  29891  63252D-04 
8.61367  05404  57314D-03 
6.73338  50887  03690D-02 
2.51448  87747  33840D-01 
6.34184  55737  37690D-01 
1.24840  40550  83 152D+00 
2.06568  80319  53401D+00 
3.00919  93586  62542D+00 
4.00041  62696  90208D+00 
5.00000  00000  OOOOOD+OO 

1.01695  09859  48944D-03 
2.29467  06865  17670D-02 
1.07665  79680  22888D-01 
2.73457  76624  65576D-01 
4.97881  55919  24992D-01 
7.25620  89195  65360D-01 
8.95263  86903  20078D-01 
9.77815  74653  81624D-01 
9.98339  07813  99277D-01 
9.99991  63424  08948D-01 
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Table  7.  (continued) 


O a Vi 

12. 0 8 5.71021  84272  06990D-04 

1.54042  43511  15548D-02 
8.83424  84071  96555D-02 
2.82446  20545  09770D-01 
6.57486  98923  05580D-01 
1.24654  10609  77993D+00 
2.03921  84951  30811D+00 
2.97933  34870  49800D+00 
3.98577  25953  93049D+00 
4.99724  08043  11428D+00 

5.99986  87939  5 1190D+00 

7.00000  00000  00000D+00 

14.0  9 3.41982  14602  49725D-04 

9.29659  34301  87960D-03 
5.40621  47717  55252D-02 
1.76394  50965  08648D-01 
4.21848  66056  53738D-01 
8.27402  28958  84040D-01 
1.41028  75856  37014D+00 
2.16099  75052  38 153D+00 
3.04350  47493  58223D+00 
4.00569  25790  69439D+00 
4.99973  27079  05968D+00 

5.99987  51919  71098D+00 
6.99999  45605  68667D+00 

8.00000  00000  OOOOOD+OO 

16.0  10  2.15843  89882  80793D-04 

5.89843  27437  09196D-03 
3.46279  59568  96131D-02 
1.14558  64950  70213D-01 
2.79034  42188  56415D-01 
5.60011  37986  53321D-01 
9.81409  12428  83119D-01 
1.55359  48539  74655D+00 
2.27017  91140  36658D+00 
3.10823  46017  15371D+00 
4.03293  08939  96553D+00 
5.00680  32702  28157D+00 
6.00081  54667  35179D+00 
7.00004  50350  79542D+00 

8.00000  07389  23901D4-00 

9.00000  00000  00000D+00 


Ui 

2.92101  89269  12141D-03 
3.43113  06112  56885D-02 
1.22466  94956  38615D-01 
2.76110  82420  22520D-01 
4.79780  96430  10337D-01 
6.96655  56772  71379D-01 
8.79007  79419  72658D-01 
9.86862  24492  94327D-01 
1.01514  23896  88201D+00 
1.00620  97126  32210D+00 
1.00052  88299  22287D+00 
1.00000  23977  96838D+00 
1.75095  72432  02047D-03 
2.08072  65842  87380D-02 
7.58683  06164  33430D-02 
1.76602  05266  71851D-01 
3.20662  43620  72232D-01 
4.93440  52905  538 12D-01 
6.70749  70306  98472D-01 
8.24495  90253  66557D-01 
9.31464  67421  62802D-01 
9.84576  84431  63154D-01 
9.99285  27691  54770D-01 
1.00027  31129  57723D+00 
1.00002  28574  02321D+00 
1 .00000  00814  05 1 80D+00 
1.10580  48735  01181D-03 
1.32449  99447  07956D-02 
4.89984  23075  92144D-02 
1.16532  61928  68815D-01 
2.17858  66931  94957D-01 
3.48176  60169  45031D-01 
4.96402  79159  11545D-01 
6.46902  61896  23831D-01 
7.82368  89717  83889D-01 
8.87777  24458  93361D-01 
9.55166  50770  35583D-01 
9.87628  55797  41800D-01 
9.97992  91838  63017D-01 
9.99847  06206  34641D-01 
9.99996  28916  45340D-01 
9.99999  99468  93169D-01 
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Table  8.  The  nodes  v\,  ...  ,vj  and  weights  u\, ...  ,uj  for  the  quadrature  rule 

SJnkah(g ) = h £'=1  Ui  g(Vih)  + h E-=o  S(ah  + ih)  + h EL i wi  g(l  - xth),  with 
h = (n+a+b—  l)_1,forg(;t)  = </>(x)  log  x + \J/(x),  with  0 and  \j/  regular  functions. 
The  error  is  of  order  0(hl  log  h).  The  nodes  x\, . . . ,x^  and  weights  wi, ...  ,Wk  are 
found  in  Table  6. 


/ 

a 

Vi 

Ui 

2 

1 

1.59154  94309  18953D-01 

5.00000  00000  00000D— 01 

3 

2 

1.15039  58119  72836D-01 
9.36546  45279  49632D-01 

3.91337  37887  53340D-01 
1.10866  26211  24666D+00 

4 

2 

2.37964  72841  18974D-02 
2.93537  07415  01914D-01 
1.02371  51242  51890D+00 

8.79594  26755  93887D-02 
4.98901  71529  13699D-01 
9.13138  85795  26912D-01 

5 

3 

2.33901  30272  03800D-02 
2.85476  49313  11984D-01 
1.00540  33272  20700D+00 
1.99497  03039  94294D+00 

8.60973  65561  58105D-02 
4.84701  96854  17959D-01 
9.15298  88691  23725D-01 
1.01390  17789  84250D+00 

6 

3 

4.00488  41949  26570D-03 
7.74565  53733  36686D-02 
3.97284  99935  23248D-01 
1.07567  33529  15104D+00 
2.00379  69271  11872D+00 

1.67187  96911  47102D-02 
1.63695  83714  47360D-01 
4.98185  65697  70637D-01 
8.37226  62455  78912D-01 
9.84173  08440  88381D-01 

8 

5 

6.53181  57085  67918D-03 
9.08674  45846  57729D-02 
3.96796  65333  75878D-01 
1.02785  66405  25646D+00 
1.94528  85929  09266D+00 
2.98014  79338  89640D+00 
3.99886  13499  51123D+00 

2.46219  41989  95203D-02 
1.70131  58668  54178D-01 
4.60925  63586  50077D-01 
7.94729  11486  21894D— 01 
1.00871  04143  37933D+00 
1.03609  36497  26216D+00 
1.00478  76565  33285D+00 

10 

6 

1.17508  93812  27308D-03 
1.87703  41298  31289D-02 
9.68646  83914  26860D-02 
3.00481  86680  02884D-01 
6.90133  15571  73356D-01 
1.29369  57380  83659D+00 
2.09018  77297  98780D+00 
3.01671  93131  49212D+00 
4.00136  97478  72486D+00 
5.00002  56617  93423D+00 

4.56074  68820  84207D-03 
3.81060  63223  84757D-02 
1.29386  49972  89512D-01 
2.88436  03814  08835D-01 
4.95811  19143  44961D-01 
7.07715  46005  94529D-01 
8.74192  43652  85083D-01 
9.66136  19865  15218D-01 
9.95788  78660  78700D-01 
9.99866  57874  23845D-01 
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Table  8.  (continued) 


l a Vi 

Y2  7 1.67422  36826  68368D-03 

2.44111  00950  09738D— 02 
1.15385  12974  29517D-01 
3.34589  84904  80388D-01 
7.32974  05318  07683D-01 
1.33230  50485  25433D+00 
2.11435  87523  25948D+00 
3.02608  45496  553 18D+00 
4.003 1 6 630 1 2 92590D+00 
5.00014  11700  55870D+00 
6.00000  10024  41 859D+00 
14  9 9.30518  23685  45380D-04 

1.37383  24584  34617D-02 
6.63075  27607  79359D-02 
1.97997  13976  22003D-01 
4.50431  35038  16532D-01 
8.57188  86311  01634D-01 
1.43450  52296  17112D+00 
2.17517  78341  37754D+00 
3.04795  50683  86372D+00 
4.00497  49068  13428D+00 
4.99852  59018  20967D+00 
5.99952  30151  16678D+00 
6.99996  36178  83990D+00 

7.99999  94881  30134D+00 

16  10  8.37152  98320  141 13D-04 

1.23938  27255  42637D-02 
6.00929  07857  39468D-02 
1.80599  12496  01928D-01 
4.14283  25990  2803 ID-01 
7.96474  77311  12430D-01 
1.34899  38824  67059D+00 
2.07347  16602  64395D+00 
2.94790  49390  31494D+00 
3.92812  92522  48612D+00 
4.95720  30865  63112D+00 
5.98636  01139  77494D+00 
6.99795  77047  91519D+00 
7.99988  87575  24622D+00 

8.99999  87543  06120D+00 


«£ 

6.36419  07807  20557D-03 

4.72396  41432  87529D-02 
1.45089  11583  85963D-01 
3.02165  94707  85897D-01 
4.98427  07397  15340D-01 
6.97121  37951  76096D-01 
8.57729  56227  57315D-01 
9.54413  65543  51155D-01 
9.91993  80527  76484D-01 
9.99462  18758  22987D-01 
9.99993  44080  92805D-01 
3.54506  06447  80164D-03 
2.68151  40315  76498D-02 
8.50409  20350  93420D-02 
1.85452  62166  43691D-01 
3.25172  43748  83192D-01 
4.91155  37472  60108D-01 
6.62293  34173  69036D-01 
8.13725  45788  405 10D-01 
9.23559  55149  44 174D-01 
9.82160  99237  44658D-01 
1.00004  73945  96121D+00 
1 .00090  93366  93954D+00 
1.00011  95342  83784D+00 
1.00000  28357  46089D+00 
3.19091  90866  26234D-03 
2.42362  13804  26338D-02 
7.74013  55216  53088D-02 
1.70488  94202  86369D-01 
3.02912  34785  11309D-01 
4.65222  08349  14617D-01 
6.40148  96370  96768D-01 
8.05121  29461  81061D-01 
9.36241  19456  98647D-01 
1.01435  97753  69075D+00 
1.03516  77210  53657D+00 
1.02030  86249  84610D+00 
1.00479  83974  41514D+00 
1.00039  50173  52309D+00 
1.00000  71494  22537D+00 


